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Abstract. This paper describes U2DE, a finite-volume code 
that numerically solves the Euler equations. The code was 
used to perform multi-dimensional simulations of the gradual 
opening of a primary diaphragm in a shock tube. Erom the 
simulations, the speed of the developing shock wave was 
recorded and compared with other estimates. The ability of 
U2DE to compute shock speed was confirmed by comparing 
numerical results with the analytic solution for an ideal shock 
tube. 

Eor high initial pressure ratios across the diaphragm, 
previous experiments have shown that the measured shock 
speed can exceed the shock speed predicted by one-dimen- 
sional models. The shock speeds computed with the present 
multi-dimensional simulation were higher than those esti- 
mated by previous one-dimensional models and, thus, were 
closer to the experimental measurements. This indicates that 
multi-dimensional flow effects were partly responsible for 
the relatively high shock speeds measured in the experi- 
ments. 

Key words: Unstructured-grids, Solution-adaptive remesh- 
ing, Numerical dissipation. Shock speeds. Shock tubes flow 



Nomenclature, Units 

A: cell area in the (x,y)-plane, m^ 
a: temperature dependent function within Redlich- 
Kwong equation of state, m^/kg • s 
oq: constant within Redlich-Kwong equation of state, 
m^/kg • s 

b: constant within Redlich-Kwong equation of state, 

mVkg 

CS: contact surface 

Cv : specific heat at constant volume, J/kg-K 
E: total energy (internal + kinetic), J/kg 
e: specific internal energy, J/kg 
F; array of flux terms 

Correspondence to: P.A. Jacobs 



n; unit normal vector 

P: pressure. Pa 

Q; array of source terms 

R\ specific gas constant, J/kg-K 

S\ surface 

T: temperature, K 

Tc'. critical temperature, K 

Tp\ reduced temperature 

t: time, seconds 

U: array of conserved quantities 
u; flow velocity, m/s 
u\ flow speed in x-direction, m/s 
v\ flow speed in y-direction, m/s 
a: noise filter coefficient 
p: density, kg/m^ 

7: ratio of specific heats 
a: Courant number 
d\ cell volume, m^ 

d'\ volume per radian for an axisymmetric cell, m^ 

Subscripts and superscripts 
L; left 
R; right 

1: driven-gas initial state 
4; driver-gas initial state 
n: time level 



1 Introduction 

Shock tunnels and expansion tubes are used to generate high 
energy flows for the ground testing of hypersonic vehicles. 
Elow in each facility is usually initiated by the rupture of 
a primary diaphragm that separates a high pressure driver 
gas and a low pressure driven gas. High pressure driver gas 
expands into the driven section, compressing and accelerat- 
ing the lower pressure driven gas. A shock wave develops 
within the driven gas and propagates along the shock tube. 
The strength of shock wave which is measured by the pres- 
sure ratio or shock speed determines the flow conditions 
behind the shock and subsequently in the test section. 

The shock speed in an ideal shock tube, with instan- 
taneous diaphragm removal and negligible viscous effects. 
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White’s Models 

Ikui’s models 

° Millers Jones 

U2DE- constant area shock tube 
° U2DE - Langley expansion tube 

Fig. 1. Maximum Mach number of developed shock wave in a 
shock tube versus initial pressure ratio. Experimental data by Miller 
and Jones (1975) is compared with the theories of White (1958) 
and Ikui et al. (1969). The Mach number of the shock wave has 
been normalized by the Mach number of the shock wave in an 
ideal shock tube for the same initial conditions 

can be determined by solving the unsteady one-dimensional 
Euler equations. The ideal shock speed has been reported 
by Duff (1959), to overestimate the shock speed in long, 
thin driven tubes where viscous effects are significant. Con- 
versely, measured shock speeds (White 1958; Miller and 
Jones 1975; Huber 1958; Nagamatsu et al. 1959; and Whit- 
field et al. 1966) in “high-performance” and larger diameter 
shock tubes can exceed the ideal shock speed by up to 20% 
when the initial pressure ratio across the diaphragm exceeds 
10^ as shown in Fig. 1. The higher-than-ideal shock speeds 
can be partially explained by considering the wave processes 
which occur during the gradual opening of a diaphragm. 



1.1 Previous work 

White (1958) developed a theory based on shock formation 
from compression waves. The model assumes that unsteady 
isentropic compression waves are formed in the driven gas as 
the diaphragm gradually opens. The compression waves are 
then assumed to coalesce into a shock wave at some distance 
downstream from the diaphragm. An upstream-facing expan- 
sion is formed to match the flow conditions. This model can 
predict higher maximum shock speeds than the ideal shock 
tube model, but it fails to predict the shock front acceler- 
ation which has been observed in experiments. As an im- 
provement to the model of White (1958), Ikui et al. (1969) 
developed a multi-stage model. They assumed that a series 



of compression waves produced by the gradual opening of 
the diaphragm, can be divided into a finite number of groups 
of compressions. A group of compression waves coalesce at 
the same point and the shock front generated by the first 
group is successively accelerated by the other groups. This 
model can predict slightly higher maximum shock speeds 
than the model of White (1958) as shown in Fig. 1. 

Zeitoun et al. (1979) performed a one-dimensional com- 
putation using the method of characteristics, in which the 
finite opening time of the diaphragm and boundary layer ef- 
fects were taken into account. The finite opening time of the 
diaphragm was found to induce a strong shock acceleration 
followed by a slow deceleration, and the maximum com- 
puted shock speed was close to the value predicted by the 
theory of White (1958). When the effects of the boundary 
layer were neglected, the shock still decelerated after the 
initial acceleration but its velocity remained higher than the 
ideal shock speed. However, the inclusion of boundary layer 
effects caused a monotonic decrease in the shock speed to 
the value below the ideal shock speed. 

Miller and Jones (1975) measured shock-wave veloci- 
ties in the Fangley six-inch diameter expansion tube. Air, 
argon, carbon dioxide and helium were used as test gases. 
The driver gas was always helium. The shock speed mea- 
surements were made using a microwave interferometer and 
via the response of pressure transducers positioned along the 
driven section (time of arrival gauges). The maximum shock 
speeds measured exceeded the maximum shock speeds pre- 
dicted by the one-dimensional theories of White (1958) and 
Ikui et al. (1969) at high initial pressure ratios for all test 
gases except argon. 

The experimental work of Miller and Jones (1975) was 
chosen to be the reference point, because of the high quality 
and detail of the available experimental data. Figure 1 com- 
pares the normalised maximum shock Mach number (helium 
as the test gas) with the theories of White (1958) and Ikui 
et al. (1969) in which constant area tube and ratio of spe- 
cific heats 7 = 1.667 are assumed. The maximum shock 
Mach number is normalised by the shock Mach number ex- 
pected in an ideal shock tube at the same conditions. Two 
important observations can be made from Fig. 1 ; (i) the ex- 
perimental data points are significantly higher than estimated 
values from the one-dimensional theories; and (ii) the nor- 
malized shock Mach number predicted by the theories of 
White (1958) and Ikui et al. (1969) increases with initial 
pressure ratio. 

Miller and Jones (1975) suggested that the higher shock 
speeds were caused by a combination of mechanisms in- 
cluding heating of driver gas during pressurisation, effects 
of the finite opening time, and multi-dimensional effects. 
The multi-dimensional nature of the flow resulting from a 
gradually opening diaphragm is examined here to determine 
if it contributes significantly to the higher-than-expected ex- 
perimental shock speeds. 

Multi-dimensional simulations of the gradual opening 
of a diaphragm have been performed previously (Sato- 
fuka 1970; Outa et al. 1975; Cambier et al. 1992; Vasil’ev 
and Danil’chuk 1994). These works have modelled the di- 
aphragm opening as a slit in two-dimensional flow or an 
iris in an axisymmetric geometry. The opening commenced 
in the middle and then progressed towards the shock tube 
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Fig. 3. Geometry associated with error function 



wall. Most of the work concentrated on the flow develop- 
ment and, in particular, on the structure of contact surface 
and the expansion waves in the driver gas. 

Satofuka (1970) performed a numerical study of shock 
formation in cylindrical and two-dimensional shock tubes. 
Air/Air driver-driven gas combinations were examined at 
diaphragm pressure ratios of 10, 100 and 1000. The calcu- 
lated shock speeds were similar to those of White (1958) 
and Ikui et al. (1969) at the lower initial pressure ratios. 
However, at the highest pressure ratio of 1000, a slightly 
higher shock speed (h- 0.05%) was predicted. 

Outa et al. (1975) performed experiments and two- 
dimensional simulations of a gradually opening diaphragm. 
From schlieren pictures and the numerical simulations, the 
presence of oblique shock waves interacting with a two- 
dimensional unsteady expansion was observed. It was con- 
cluded that the effects of these waves on the flow structure 
were restricted to one to two diameters downstream. The 
maximum experimentally measured shock speed within a 
100 mm square shock tube for an initial pressure ratio of 
6100 exceeded the ideal shock speed by 10% and exceeded 
the maximum shock speed predicted by the theory of Ikui 
et al. (1969) by 5%. 

Cambier et al. (1992) performed a two-dimensional ax- 
isymmetric simulation of gradual diaphragm rupture in which 
the diaphragm was modelled as an opening iris. The fol- 
lowing observations were made from the simulations; the 
primary shock becomes planar very rapidly (within two di- 
ameters from the diaphragm), a complex and unsteady flow 
structure dominated by a Mach disk is formed behind the 
contact surface (CS), and the CS itself becomes a complex 



shape. The initial shape of the CS is due to the relatively 
slow opening time of the diaphragm. The contact surface CS 
does not become planar with time, and it was suggested that 
its fate could be dominated by Rayleigh-Taylor instabilities. 

Vasil’ ev and Danil’chuk (1994) performed an inviscid 
two-dimensional simulation of shock wave formation in a 
shock tube by considering transverse diaphragm removal. 
Two main observations were made from their simulations: (i) 
jetting of the CS along the walls due to a system of oblique 
shock waves; and (ii) fragmentation of the secondary shock 
which occurred because a pocket of hot unexpanded gas at 
the wall changed the effective area of the tube. The resulting 
flow is analogous to the flow through a Laval nozzle. 



1.2 Scope of the current work 

The results from two-dimensional axisymmetric simulations 
of a gradually opening diaphragm with high initial pressure 
ratios are presented here. The simulations were performed 
using a finite-volume code, U2DE, which solves the Euler 
equations. The diaphragm opening is modelled as an iris, 
similar to the model proposed by Cambier et al. (1992). 
Particular attention is directed to the variation of the speed 
of the developing shock. 

Experimental conditions approximating those reported 
by Miller and Jones (1975) were investigated. The mod- 
els of White (1958) and Ikui et al. (1969) fail to predict the 
maximum shock speed at these conditions. It will be shown 
that the multi-dimensional nature of flow contributed to the 
higher than expected maximum shock speed. The structure 
of the flow as it developed during and after diaphragm rup- 
ture, and in particular the contact surface, are also examined. 



2 Computational model 



A cell-centred finite-volume code, U2DE was used to solve 
numerically the Euler equations. The flow domains were 
represented as unstructured meshes of triangular cells and 
solution-adaptive remeshing was used to focus computa- 
tional effort in regions where the flow-field gradients were 
high. 

The two-dimensional Euler equations can be written as. 



^ f Vd'd+ f F-MS' = 0, 
dt A Js 



where 



U = 




F = 



puu + Pi 
pvu + Pj 
. pEu + Pu j 



( 1 ) 



( 2 ) 



The two-dimensional cells are assumed to have unit depth. 
The axisymmetric form of the Euler equations can be written 
similarily as. 



where. 



( 3 ) 
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Fig. 4a, b. Initial mesh: a for ideal shock-tube problem (P4/P1 = 10) and a solution-adapted mesh; b at time t = 6 x lO”"* s 



Q = 




( 4 ) 



The volume of the cell is expressed as volume per radian 

For calorically perfect gas, the equation of state of ideal 
gas is used, 



P = il- l)pe. 



( 5 ) 



The specific internal energy is written as, 
e = E + v^). (6) 



However, the initial pressure of the helium driver gas used 
by Miller and Jones (1975) was high enough to cause devi- 
ation from ideal gas behaviour due to van der Waals forces. 
The Redlich-Kwong equation of state (Aungier 1995) can 
be used in (7) to describe more accurately the behaviour of 
helium at the flow conditions of interest. 



pRT p^a(T) 

I — bp l + bp ^ 



( 7 ) 



a(T) = aoT-«“, 



( 8 ) 



Tr = 



T 

7;’ 




( 9 ) 

( 10 ) 



where ao= 226.20 m^/kg-s^, b = 4.1648 x 10 ^ m^/kg, C„= 
31 15.6 J/kg-K, R= 2077 J/kg-K and = 5.3 K. The specific 
heat Cy, was assumed to be constant. 

The Euler equations are applied to each triangular cell 
in the discretised form, 

f 

k=\ 

and, from known initial conditions, the flow solution in each 
cell is explicitly updated in time. Each time-step can be split 
into three parts. Eirstly, the pseudo-left and -right edge flow 
states are determined at the edges that bound each triangular 
cell. Secondly, the flux array F at each edge is determined. 
Einally, the cell-averaged conserved quantities U and the 
primary flow variables for each cell are updated. 

The pseudo-left and -right edge flow states are recon- 
structed from the cell averaged data with each primary flow 
variable being treated independently. Eor example, the left 
and right edge densities (pj^„, p„„) are constructed from the 



densities at four nearby points (p^^, p^^, p^^, p^^) as shown 
in Eig. 2. If the edge is internal to the flow field, p ^2 equals 
the density at the vertex of the left cell which is opposite 
the edge, p^i equals the density at the centre of the left cell, 
Pr\ equals the density at the centre of the right cell, and p^a 
equals the density at the vertex of the right cell opposite to 
the edge. A primary flow variable at a vertex is determined 
by summing the primary flow variable of the surrounding 
cells multiplied by a weight, and then dividing by the sum 
of the weights. The weight is equal to the inverse of the 
distance from the vertex to the centre of the cell (Batina 
1993). 

If the edge is external, the cell associated with the edge is 
defined to be the right cell. The density, pressure and internal 
energy of the left and far-left pre-interpolation values are set 
to the right and far-right cell values respectively. The left 
and far-left velocities are set to the reflected velocities of 
the right and far-right cell respectively with the edge acting 
as a mirror. 

To model the gradual opening of the diaphragm, U2DE 
has the ability to blank out (ignore) parts of the domain. The 
use of unstructured meshes made the implementation of this 
feature easy. Edges between an ignored cell and a flow cell 
are treated as a wall. 

A generalized MUSCL interpolation scheme (Anderson 
et al. 1985) is used to construct the left and right flow states 
from the pre-interpolation flow states as 

Plq = Pl\ + 

Pm = Pm + (12) 

where 

= ^[(1 - «)MM{(pli - Pl2),/3(pri - PLl)} 

H-(l H- K)MM{f3(pLl - PL 2 ),(Pm - PLl)}], 
'^)MM{(pri — pli), /3(pr2 — Pri)} 

+(1 - k)MM{P(pri - pLi), ipR2 - PRl)}]. (13) 

The minmod (MM) limiter function returns the argument 
with the minimum magnitude if both arguments have the 
same sign and returns zero otherwise. The parameter «; = 1/3 
was used giving an upwind-biased third-order interpolation 
scheme, and a compression parameter /3 = 2 was used. 

The equilibrium flux method (EEM) originally developed 
by Pullin (1979) is used to calculate the flux array from 
the left and right edge flow states. EEM is derived from 
the kinetic theory of gases, and it has been demonstrated 
by Macrossan (1989) that EEM solves the Euler equations 
with added pseudo dissipation and, in the hypersonic limit. 
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Fig. 5a-c. Comparison of numerical density profiles for an ideal 
shock-tube problem (P 4 /P 1 = 10) with analytical solution at t = 
600/is: a first order solution; b higher order solution; and c higher 
order solution with solution-adaptive remeshing 



EFM becomes an upwind scheme. The EFM flux calculation 
assumes perfect gas. For a non-perfect gas, an effective ratio 
of specific heats (14) calculated from the pseudo left and 
right edge flow states is used (Edwards 1988). 



^ 



(14) 

(15) 



The increment in flow time for each time-step is determined 
during the flux calculations and is equal to, 



At = a X minimum 



/ local wave speed \ 
\ smallest cell median J 



(16) 



where a, is the Courant number and is usually set to 0.5. 

The cell averaged conserved quantities are advanced 
from time level n to time level n + I using the predictor- 
corrector scheme 



dt 



u6) = H- 



= Z\f 



dt 



U(n+1) = u(i) + _(/\u<2) _ Z\u6)), 



(17) 



where the superscripts 6) and indicate intermediate re- 
sults. If a first order scheme is desired, only the first stage 
is used and = u(i)^ 



2.7 Solution-adaptive remeshing 

Solution-adaptive remeshing concentrates the computational 
effort at regions of interest within the flow domain. This al- 
lows better resolution of discontinuities such as shock waves 
and slip lines than would be possible with fixed-grid simu- 
lations at the same (or similar) computational expense. The 
resolution of the mesh is increased by introducing nodes to 
the mesh thereby increasing the number of cells in that re- 
gion. The resolution of the mesh can be reduced in regions 
where the solution has become smooth by removing previ- 
ously inserted nodes. 

The remeshing process comprises three stages: firstly, 
the “error indicator” is calculated for each cell and cells 
are marked for deletion, refinement or no action; the second 
step is the deletion of vertices surrounded by cells which 
have been marked for deletion; and finally, cells marked for 
refinement are split. The frequency of remeshing depends on 
the Courant number and the number of “protective layers” 
provided during refinement. A protective layer is formed by 
refining the cells adjacent to the cells marked for refinement. 
For a Courant number a = 0.5 and using one protective 
layer, remeshing was performed every five time-steps. 

The primary flow variable used to calculate the error 
indicator is density. The error indicator associated with each 
cell is determined by 

error indicator = | 1 26 — Oj — | 

i,3,k 

{^(|6- a*| H- |ci - b\) + a’^(ai + 2b + Ci)Y (18) 

iyj,k itjik 

The geometry associated with this equation is shown in 
Fig. 3 where Ui is the density at the centre of an adjacent cell, 
b is the density at the centre of the cell, and Ci is the density 
at the vertex opposite to the adjacent cell. This indicator is 
based on similar error functions developed by Fohner (1987) 
and Probert et al. (1991). 

If the error indicator is greater than 0.3 the cell is marked 
for refinement, and if less than 0.1 the cell is marked for 
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deletion. If the volume of a cell is lower than the specified 
minimum volume or area for axisymmetric cases, it is not 
marked for refinement. 

The bisection method (Rivara 1984) is used to refine the 
triangular cells. Cell refinement is achieved by inserting a 
new vertex at the midpoint of an edge and splitting the cells 
adjoining the edge. The edge chosen to be split must have 
the greatest length of all the edges of the two cells adjoining 
the edge. A refinement level is assigned to a vertex when it is 
inserted. This number is equal to the highest refinement level 
of the surrounding vertices plus one. Initially the refinement 
level of all vertices is zero. 

Only vertices associated with four or two cells (boundary 
vertex) are considered for deletion. All the cells connected 
to the vertex must be marked for deletion and the refinement 
level of the vertex must be higher than all vertices connected 
to the vertex being considered for deletion. When a vertex is 
inserted, the index numbers of the vertices of the split edge 
are stored within the data structure of the new vertex. This 
information is used to ensure that deletion is the reverse of 
a previous insertion. The retention of the local mesh allows 
for further vertex removal. 

The cell refinement for the axisymmetric simulations was 
carried out in a manner so that cell aspect ratio was main- 
tained throughout the simulation. “Numerically induced jet- 
ting” along the axis has been previously observed in ax- 
isymmetric simulations by Cambier et al. (1992), and by the 
authors. Stretching the cells in the radial direction can alle- 
viate this problem. A cell aspect ratio of three was used for 
all the axisymmetric simulations presented here. 



3 Code validation - ideal shock tube 

The accuracy of U2DE was validated by comparing nu- 
merical solutions for an idealised shock tube to the one- 
dimensional analytical solution. 

Firstly, we examine a shock tube with a low initial pres- 
sure ratio across the diaphragm (Hirsch 1990). The gas is 
assumed to be calorically perfect with 7 = 1.4. The initial 
state at time f = 0 is 

a; < 0.5 m : p 4 = 1.0 kg/m^, 

P 4 = 10^ Pa, U 4 = V 4 = 0, 
x > 0.5 m : Pi = 0.125 kg/m^, 

Pi = 10"^ Pa, ui = vi= 0. (19) 

Figure 4 shows the initial mesh and a solution-adapted mesh 
at a later time. Profiles of density from first-order and higher- 
order solutions generated without solution adaptive remesh- 
ing and a higher-order solution generated using solution 
adaptive remeshing are shown in Fig. 5. The volume of each 
cell within the initial mesh is 5.0 x 10“^ m^. The minimum 
cell volume was set to 1.0 x 10 ”^m^ and the coefficient of 
the noise filter, a in the error indicator (18) was set to 0 . 01 . 
The comparisons of the numerical solutions with the analyt- 
ical solution demonstrate that the current algorithm provides 
a satisfactory solution and that the implementation of the 
higher order interpolation and solution adaptive remeshing 
does improve the accuracy of the solution. 



We now examine the ability of U2DF to produce accu- 
rate numerical solutions when the initial pressure ratio across 
the diaphragm is high. The initial state at time f = 0 is set 
to 

a; < 5.0 m : P 4 = 35 X 10® Pa, 

T 4 = 342 K, U 4 = 0, 

X > 5.0 m : Pi = 3450 Pa, 

Ti = 297.6 K, ui = 0. (20) 

This condition approximates one of the experimental condi- 
tions used by Miller and Jones (1975), with P 4 /P 1 = 10, 145, 
and provides a harsher test of the code. Helium is used as 
driver gas as well as test gas. The high pressure of the driver 
gas (x < 5.0 m) leads to significant deviation from the 
perfect-gas model and so the Redlich-Kwong equation of 
state is used. 

Two-dimensional planar simulations were performed. 
The volume of each cell within the initial mesh was 2.0 x 
10“"^ m^. The minimum cell volume for the simulation was 
set to 1.0 X 10“® m^ and a = 0.01. Figure 6 compares 
the higher order numerical density profile with the one- 
dimensional analytical solution. The agreement is satisfac- 
tory but the shock position does appear to be incorrectly es- 
timated. Thus, the accuracy of U2DE to compute the shock 
speed was examined in more detail. 

The shock speed was computed by recording the po- 
sition of the shock wave every ten time-steps. This data 
was smoothed and then differentiated. Speeds of developing 
shocks were computed at various initial pressure ratios ( 10 , 
100 , 1 , 000 , and 10 , 000 ) across the diaphragm and compared 
to the corresponding analytical shock speeds in Fig. 7. The 
density and pressure on the right side of the tube was set 
to unity. The initial temperature of the driver and driven 
gas was the same for all cases and the gas was assumed 
to be calorically perfect with 7 = 1.667. The simulations 
used the same domain and initial mesh as the first shock 
tube problem as seen in Fig. 4. The minimum cell volume 
was set to 1.0 x 10~^ m^ and a = 0.03. At the lower initial 
pressure ratios the agreement between the computed shock 
speed and correct (analytical) shock speed is good. At higher 
initial pressure ratios, there is an initial overestimation, but 
the computed shock speed does decay to the correct shock 
speed. The magnitude of the initial overestimation and the 
distance required to decay to the analytical value increases 
with initial pressure ratio. 

It was speculated that the initial overestimation of shock 
speed when the initial pressure is high is due to numerical 
diffusion, particularly at the contact surface. A test case was 
designed such that the initial pressure ratio across the di- 
aphragm was high ( 10 '^), and the densities on either side of 
the contact surface were equal. The domain and initial mesh 
for the first shock tube problem as seen in Fig. 4 were used. 
The initial condition was, 

a; < 0.5 m : p 4 = 11.9969 kg/m^, 

P 4 = 10^ Pa, 

a; > 0.5 m : Pi = 1.0 kg/m^. 

Pi = 1.0 Pa. (21) 
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Fig. 6. Density profile for an ideal shock tube {P^,/ P{ = 10, 145) 
at f = 500/rs. Comparison of numerical solution with analytical 




X, (m) 



Fig. 7. Shock Mach number versus distance from diaphragm for an 
ideal shock tube. Numerical values at various initial pressure ratios 
are compared with analytical values. Note the diaphragm location 
is at a; = 0.5 m 



Note that this initial condition is extreme due to high initial 
temperature ratio and the large shock Mach number in Table 
1 . 

Figure 8 compares the computed density profile with the 
analytical solution at t= 8.7 ms. The position of the shock 
wave agrees with the analytical position and the distance 
required for the shock speed to settle to the analytical speed 
is significantly reduced as shown in Fig. 9. 

The computed shock speed versus distance from the di- 
aphragm for the experimental condition stated above was 
compared to the analytical shock speed as seen in Fig. 10. 
The speed of the developing shock was computed at differ- 
ent mesh resolutions (that is, minimum cell volume). The 
distance required for the shock speed to decay to the analyt- 
ical shock speed decreased with increasing mesh resolution 
(lower minimum cell volume). This is consistent with the 
reduction in numerical diffusion that is expected with in- 
creased mesh resolution. The solution adaptive remeshing 



procedure is therefore an important component of the code 
when trying to estimate accurately the shock speeds of a 
high performance shock tube. The overestimation of the pri- 
mary shock speed by the fixed-grid Navier-Stokes code can 
be seen in a previous study of the NASA Langley expansion 
tube by Jacobs (1994) as seen in Fig. 7. 

4 Simulation of gradual diaphragm opening 

Numerical simulations of flow through a gradually open- 
ing primary diaphragm are now presented. The geometry 
of the domain is based on the NASA Langley expansion 
tube (Miller and Jones 1975) operated with the primary di- 
aphragm only. The diameter of the driver section was 165.1 
mm and its length was 2.44 m. The diameter of the driven 
tube was 152.4 mm. The transition from driver tube diam- 
eter to driven tube diameter occurred after the diaphragm 
location and extended over a length of 190.5 mm. Although 
the diaphragm section was square in cross-section and the 
transitional area change piece went from square to circular, 
the geometry for the current simulations was assumed to be 
axisymmetric. 

Three experimental conditions used by Miller and Jones 
(1975) were examined in Table 2. The initial driver gas state, 

X < 2.44 m : P 4 = 35 X 10® Pa, 

T 4 = 342 K, U 4 = V 4 = 0, (22) 

was the same for all conditions. Both driver and driven gases 
were helium and the Redlich-Kwong equation of state was 
used to describe the gas behavior. The diaphragm opening 
was modelled as a dilating iris with an opening time of 
200 /is. The minimum cell area, in the {x, j/)-plane, of the 
initial mesh was 9.26 x 10““* m^, except at the diaphragm 
section. The diaphragm was created by refining the cells 
at the diaphragm location until the cell area was less than 
5.0 X 10^^ m^. A thin strip of cells, 4.24 mm thick was 
chosen to represent the diaphragm. These cells were initially 
ignored. During the simulation, the status of the ignored cells 
was changed to flow cells so that the flow area open was 
proportional to the elapse of time from its opening. The 
minimum cell area for the simulation was set to 5.0 x 10“^ 
m^ and a = 0 . 02 . 



4.7 Flow development 

The time history of the density and pressure contours for 
the Langley expansion tube (P 4 /P 1 = 10,145) is shown in 
Figs. 11-13. The initial shape of the shock wave is spher- 
ical until it reflects at the tube wall. The shock front and 
the reflected (transverse) waves interact causing the shock 
front to become planar within a relatively short distance 
and the transverse waves become weaker. These observa- 
tions of the shock front formation are similar to those made 
from previous numerical (Cambier et al. 1992; Vasil’ev and 
Danil’chuk 1994) and experimental work (Henshall 1957; 
cited by Miller and Jones 1975). 

Because the initial opening of the diaphragm occurs at 
the centre of the diaphragm, the initial shape of the contact 
surface is convex when viewed from the downstream end 
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Table 1. Summary of test cases at various initial pressure ratios where the temperatures on 
either side of the diaphragm are equal unless specified otherwise 



Pressure 

ratio 


Shock 

Mach number 


Contact surface 
density ratio 


Flow 
time (s) 


CPU 

(s) 


Final number 
of cells 


Time 

steps 


10 


1.5520 


2.5942 


0.20 


702 


2,029 


3,393 


100 


2.1945 


7.3295 


0.14 


3,685 


7,552 


3,490 


1,000 


2.7844 


21.153 


0.11 


6,615 


11,089 


3,517 


10,000 


3.2491 


59.473 


0.095 


9,565 


16,515 


3,640 


10,000“ 


35.611 


1.0000 


0.0087 


788 


1,954 


3,539 



“ Initial temperature ratio is 833.55 




Fig. 8. Ideal shock tube with a high initial pressure ratio (lO”*), but 
with no contact surface discontinuity. Comparison of numerical 
solution (o) with analytical solution (dashed line) at t= 8.7 ms. 
Minimum cell volume of 1.0 x 10~’ m^ and a = 0.01 




Fig. 9. Computed shock Mach number (solid line) versus distance 
compared with analytical Mach number (dashed line) for an ideal 
shock tube with high initial pressure ratio of 10,000 but no contact 
surface discontinuity. Note the diaphragm location is at a; = 0.5 m 




Fig. 10. Shock speed versus distance from diaphragm for an ideal 
shock tube. The initial condition is the same as an experimental 
condition of Miller and Jones (1975). Computed speeds at vari- 
ous mimimum cell volumes (m^) are compared with the analytical 
speed 



of the tube. At approximately 100 /is the interaction of the 
radially expanding driver gas with the tube wall causes an 
oblique upstream-facing shock to develop. This shock redi- 
rects the flow along the wall and causes the density and the 
pressure of the flow behind the contact surface to be higher 
at the wall than at the central part. This region of higher 
pressure gas accelerates the contact surface at the wall rel- 
ative to the centre. The contact surface eventually becomes 
concave when viewed from the downstream end. The evolu- 
tion of the contact surface is similar to that observed in the 
previous numerical studies. Note that the study of Cambier 
et al. (1992) took into account viscous effects which slowed 
down the contact surface at the walls, but the jetting of the 
contact surface near the walls relative to the center of the 
tube was evident. 

Cambier et al. (1992) discussed the effects that Rayleigh- 
Taylor instabilities may have on the development of the con- 
tact surface. Taylor (1950) showed that a contact surface be- 
tween two fluids, experiencing an acceleration perpendicular 
to their interface is stable if the heavier fluid is pushing the 
lighter fluid and unstable if the opposite is true. Since, for the 
present conditions, the density of the expanded driver-gas, is 
greater than the shock-processed driven gas, the contact sur- 
face will be stable during an acceleration phase and unstable 
during a deceleration phase. 
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Fig. 11. Time history at 20 - 80 fis of density (lop) and pres- 
sure (bottom) contours for the numerical simulation of the primary 
diaphragm opening within the NASA Langley facility 



4.2 Shock speed 

The computed shock speeds as functions of distance down- 
stream from the diaphragm for the initial conditions stated 
in Table 2 are compared with the experimentally measured 
shock speeds shown in Fig. 14. The maximum experimen- 
tal shock speeds exceed the computed speeds for all cases, 
however, the experimental and computed profiles are simi- 
lar in that both exhibit an acceleration phase followed by a 
deceleration phase. Note that the computed profile has a de- 




X, (m) 

Fig. 12. Time history at 100 - 160 ps of density (top) and pres- 
sure (bottom) contours for the numerical simulation of the primary 
diaphragm opening within the NASA Langley facility 



celeration phase even though viscous effects are not included 
which was also noted by Zeitoun et al. (1979). 

The grid convergence of the computed shock speeds was 
examined in Fig. 15, and only occurred for P 4 /P 1 = 1014.5. 
This is similar to results for the ideal shock tube as dis- 
cussed in Section 3, where the shock speed converged to 
the ideal value for P 4 /P 1 < 1 , 000. The simulation with 
the highest mesh resolution for P 4 /P 1 = 10, 145 required 22 
days of computation time (on a SGI R8,000 processor; 85 /is 
per cell per predictor-corrector time-step) and a higher res- 
olution simulation could not be obtained with the available 
computing resources. 

Due to the uncertainty of the computed maximum shock 
speed when the initial pressure ratio was high, axisymmetric 
simulations of gradual diaphragm opening were performed 
at lower initial pressure (10, 100, 1,000) ratios. The simu- 
lations were for a constant diameter tube (152.4 mm) with 
a diaphragm opening time of 200 /is. The gas was assumed 
to be perfect helium (7 = 1.667). The driven gas fill condi- 
tion was Pi = 1.0 X 10“^ kg/m^, Pi = 623.1 Pa. The driver 
gas and driven gas were assumed to have the same initial 
temperature (300 K). 

The computed speeds of the shock waves versus dis- 
tance downstream from the diaphragm location are shown 
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Fig. 13. Time history at 180 - 240 ps of density {top) and pres- 
sure (bottom) contours for the numerical simulation of the primary 
diaphragm opening within the NASA Langley facility 

in Fig. 16. Grid convergence was only achieved for Pa,jP\ = 
1000. For P 4 /P 1 = 10 and 100 the difference between the 
highest and middle resolution is greater than the difference 
between the middle and lowest resolution. However, the 
changes are small for P\jP\= 10. The reason for this is 
presently unknown, but we suspect numerical jetting, which 
has been demonstrated to become worse for higher resolu- 
tion axisymmetric simulations (Cambier et al. 1992). 

The maximum shock speed can occur within half a me- 
tre downstream of the diaphragm location as seen in Fig. 16. 
This is due to the interaction of the spherical shock wave 
with the shock tube wall. The maximum shock speed re- 
ferred to by this paper is the maximum developed shock 
speed after the shock has become planar and initial tran- 
sients have settled. 

The maximum shock speeds for the simulations of grad- 
ual diaphragm opening are compared to theoretical and ex- 
perimental shock speeds as seen in Fig. 1. Note that grid 
convergence was only achieved for Pa/P\ = 1,000 and 
1014.5. Despite this, the trend is consistent; the shock speed 
in the axisymmetric shock tube with a gradually opening 
diaphragm is greater than the speed predicted by various 
one-dimensional theories (White 1958; Ikui et al. 1969) for 
the same initial conditions. It is believed that this is related to 
the oblique upstream facing shock that temporarily appears 
downstream of the diaphragm which raises the entropy of 
the driver gas. Zeitoun et al. (1979) showed, using a one- 
dimensional model that, if an upstream facing normal shock 
exists downstream of the expansion, the speed of the shock 
wave can transiently exceed the ideal value. The theories of 
White (1958) and Ikui et al. (1969) do not consider this up- 
stream facing shock. The idea of increasing the entropy of 
the driver gas to generate faster shocks has been studied by 
Bogdanoff (1990) and Kendall et al. (1996), and it appears 
that similar entropy raising mechanisms are operating here. 




a X, (m) 




b X, (m) 




c X, (m) 

Fig. 14a-c. Computed shock speed versus distance from diaphragm 
location within NASA Langley expansion tube assuming gradual 
diaphragm opening. Experimental data (Miller and Jones 1975) 
and the ideal values are included. The initial pressure ratios are: a 
1014.5; b 10,145; and c 101,450 
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a X, (m) 





Fig. 15a-c. Computed shock speeds within NASA Langley expan- 
sion tube assuming gradual diaphragm opening. The shock speeds 
were computed from simulations at various minimum cell areas 
(m^). The initial pressure ratios are: a 1014.5; b 10,145; and c 
101,450 



The computed shock speeds obtained via the multi- 
dimensional model, although higher than the one-dimension- 
al shock speeds, are less than the experimental values of 
Miller and Jones (1975). There are a number of possible 
reasons for this. The simulations did not include the viscous 
and turbulent mixing that occurs at the contact surface. The 
temperature of the expanded gas can be very low (16 K for 






Fig. 16a-c. Computed shock speeds in a constant area shock tube 
with gradual diaphragm opening: a P 4 /P 1 = 10; b P 4 /P 1 = 100; 
and c P 4 /P 1 = 1 , 000. The Mach number of the shock wave has 
been normalized by the Mach number of the shock wave for an 
ideal shock tube with the same initial conditions. The shock speeds 
were computed from simulations at various minimum cell areas 
(m") 
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Table 2. Pressure and temperature of driven gas and maximum 
shock speeds from experiments by Miller and Jones (1975). The 
driver conditions were P = 35 MPa and T = 342 K 



Pressure (kPa) 


Temp (K) 


Max shock speed (m/s) 


34.5 


297.0 


3,490 


3.45 


297.6 


4,206 


0.345 


297.6 


4,511 



Pa,/P\= 1,000), and at these temperatures, the behaviour of 
the gas cannot be considered ideal. Also the opening of the 
primary diaphragm via petaling is fully three-dimensional 
and the current simulation is not modelling this process fully. 

As a final note, it has been shown that numerical diffu- 
sion does affect the computed shock speed for shock tube 
simulations when Pa/P\ > 1, 000. Considering this, the dif- 
fusive EFM flux calculator may appear to be a poor choice. 
However, an approximate Riemann solver (Jacobs 1992) 
which is less dissipative was also tried. When the initial 
pressure was high, this less dissipative method generated 
unacceptable levels of noise in the solution, particularly be- 
hind the shock. This phenomena, may be related to odd-even 
decoupling as discussed by Quirk (1994). 

5 Conclusion 

The finite-volume code U2DE solves the Euler equations for 
compressible flow and can be used to model shock tube flow 
and accurately compute shock speeds when the initial pres- 
sure ratio is low. When the initial pressure ratio is high the 
flow is difficult to resolve because of the large density ratio 
at the contact surface where significant numerical diffusion 
occurs. However, solution-adaptive remeshing can be used 
to control the error and obtain reasonable estimates for the 
shock speed. 

Axisymmetric simulations of the flow through a slowly- 
opening diaphragm were performed. The structure of the 
developing flow was similar to flows observed by previous 
experimental and numerical work, and the maximum speeds 
of the primary shock wave for the multi-dimensional simu- 
lations of diaphragm opening exceeded the speeds predicted 
by previous one-dimensional theories (White 1958; Ikui et 
al. 1969). It is possible that one of the mechanisms behind 
the increase in the shock speed is an entropy rise through 
the oblique shock structure which exists temporarily down- 
stream of the diaphragm while it is opening. The mechanism 
is a multi-dimensional flow effect that can be captured only 
in two- or three-dimensional simulations. 

Acknowledgements. The computer simulations were run on the 
University of Queensland High Performance Computing Facility. 
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